#### FIGURE 1: POOLED TREATMENT EFFECTS
#### Meta-analysis of economic and humanitarian narrative effects

rm(list = ls())
source("./2_code/00_setup.R")

#### STUDY 1 ####

data <- fread(paste0(data_path, "data_study1.csv"), header = TRUE)
data$worked <- as.factor(data$worked)

# Narratives on policy
lin_entry_economic <- lm_lin(exclusion_s ~ treat_economic, covariates = ~ 
                               edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_entry_humanitarian <- lm_lin(exclusion_s ~ treat_humanitarian, covariates = ~  
                                   edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])

# Narratives on prejudice
lin_prejudice_economic <- lm_lin(prejudice_s ~ treat_economic, covariates = ~ 
                                   edad + worked + male, data = data[(treat=='economic' | treat=='control')])
lin_prejudice_humanitarian <- lm_lin(prejudice_s ~ treat_humanitarian, covariates = ~  
                                       edad + worked + male, data = data[(treat=='humanitarian' | treat=='control')])


#### STUDY 2 ####

data2 <- fread(paste0(data_path, "data_study2.csv"), header = TRUE)
data2$worked <- as.factor(data2$worked)

# Support for deportation
lin_entry_economic2 <- lm_lin(exclusion_s ~ treat_economic, covariates = ~  
                                edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_entry_humanitarian2 <- lm_lin(exclusion_s ~ treat_humanitarian, covariates = ~ 
                                    edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])

# Prejudice against migrants
lin_prejudice_economic2 <- lm_lin(prejudice_s ~ treat_economic, covariates = ~ 
                                    edad + worked + male, data = data2[(treat=='economic' | treat=='control')])
lin_prejudice_humanitarian2 <- lm_lin(prejudice_s ~ treat_humanitarian, covariates = ~ 
                                        edad + worked + male, data = data2[(treat=='humanitarian' | treat=='control')])


#### STUDY 3 ####

data3 <- fread(paste0(data_path, "data_study3.csv"), header = TRUE)
data3$employed <- as.factor(data3$employed)

# Support for deportation
lin_entry_economic3 <- lm_lin(exclusion_s ~ treat_economic, covariates = ~ 
                                edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_entry_humanitarian3 <- lm_lin(exclusion_s ~ treat_humanitarian, covariates = ~ 
                                    edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])

# Prejudice against migrants
lin_prejudice_economic3 <- lm_lin(prejudice_s ~ treat_economic, covariates = ~ 
                                    edad + employed + male, data = data3[(treat=='economic' | treat=='control')])
lin_prejudice_humanitarian3 <- lm_lin(prejudice_s ~ treat_humanitarian, covariates = ~ 
                                        edad + employed + male, data = data3[(treat=='humanitarian' | treat=='control')])


#### META ANALYSIS ####

### Economic treatment, entry outcome ###

out_economic <- rbind(
  tidy(lin_entry_economic) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_entry_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_entry_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

dataset_names <- c("data", "data2", "data3")

calculate_total_observations <- function(dataset) {
  dataset %>%
    filter(treat == 'economic' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(exclusion_s)))))
}

total_observations <- lapply(dataset_names, function(name) {
  calculate_total_observations(get(name))
})

obs_eco <- unlist(total_observations)

types_economic <- data.frame(
  treatment = rep('economic', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_eco
)

out_economic <- cbind(out_economic, types_economic)

m.gen_eco <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_economic,
  common = FALSE,
  random = TRUE,
  method.tau = "REML",
  n.e = n_total,
  title = "economic_treatment_entry"
)

summary(m.gen_eco)


### Humanitarian treatment, entry outcome ###

out_humanitarian <- rbind(
  tidy(lin_entry_humanitarian) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_entry_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_entry_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

dataset_names_hum <- c("data", "data2", "data3")

calculate_total_observations <- function(dataset) {
  dataset %>%
    filter(treat == 'humanitarian' | treat == 'control') %>%
    summarise(total_observations = sum(complete.cases(across(c(treat_humanitarian)))))
}

total_observations_hum <- lapply(dataset_names_hum, function(name) {
  calculate_total_observations(get(name))
})

obs_hum <- unlist(total_observations_hum)

types_humanitarian <- data.frame(
  treatment = rep('humanitarian', 3),
  Study = c('Study 1', 'Study 2', 'Study 3'),
  n_total = obs_hum
)

out_humanitarian <- cbind(out_humanitarian, types_humanitarian)

m.gen_hum <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_humanitarian,
  common = FALSE,
  random = TRUE,
  method.tau = "REML",
  n.e = n_total,
  title = "humanitarian_treatment_entry"
)

summary(m.gen_hum)


### Economic treatment, feelings outcome ###

out_economic_feel <- rbind(
  tidy(lin_prejudice_economic) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_prejudice_economic2) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error),
  tidy(lin_prejudice_economic3) %>% filter(term=='treat_economic') %>% dplyr::select(estimate, std.error)
)

out_economic_feel <- cbind(out_economic_feel, types_economic)

m.gen_eco_feel <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_economic_feel,
  common = FALSE,
  random = TRUE,
  method.tau = "REML",
  n.e = n_total,
  title = "economic_treatment_feelings"
)

summary(m.gen_eco_feel)


### Humanitarian treatment, feelings outcome ###

out_humanitarian_feel <- rbind(
  tidy(lin_prejudice_humanitarian) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_prejudice_humanitarian2) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error),
  tidy(lin_prejudice_humanitarian3) %>% filter(term=='treat_humanitarian') %>% dplyr::select(estimate, std.error)
)

out_humanitarian_feel <- cbind(out_humanitarian_feel, types_humanitarian)

m.gen_hum_feel <- metagen(
  TE = estimate,
  seTE = std.error,
  studlab = Study,
  data = out_humanitarian_feel,
  common = FALSE,
  random = TRUE,
  n.e = n_total,
  method.tau = "REML",
  title = "humanitarian_treatment_feelings"
)

summary(m.gen_hum_feel)


#### PLOTTING ####

# Extract data for economic treatment, deportation outcome
total_observations <- sum(m.gen_eco$n.e)
forest_data_eco <- data.frame( 
  Study = c(m.gen_eco$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_eco$TE, m.gen_eco$TE.random),
  Observations = c(m.gen_eco$n.e, total_observations),
  SE = c(m.gen_eco$seTE, m.gen_eco$seTE.random),
  Lower95 = c(m.gen_eco$lower, m.gen_eco$lower.random),
  Upper95 = c(m.gen_eco$upper, m.gen_eco$upper.random),
  Lower90 = c(m.gen_eco$TE - (1.645 * m.gen_eco$seTE), m.gen_eco$TE.random - (1.645 * m.gen_eco$seTE.random)),
  Upper90 = c(m.gen_eco$TE + (1.645 * m.gen_eco$seTE), m.gen_eco$TE.random + (1.645 * m.gen_eco$seTE.random)),
  Narrative = "Economic Narrative",
  Outcome = "Deportation"
)

# Extract data for economic treatment, feelings outcome
total_observations <- sum(m.gen_eco_feel$n.e)
forest_data_eco_feel <- data.frame( 
  Study = c(m.gen_eco_feel$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_eco_feel$TE, m.gen_eco_feel$TE.random),
  Observations = c(m.gen_eco_feel$n.e, total_observations),
  SE = c(m.gen_eco_feel$seTE, m.gen_eco_feel$seTE.random),
  Lower95 = c(m.gen_eco_feel$lower, m.gen_eco_feel$lower.random),
  Upper95 = c(m.gen_eco_feel$upper, m.gen_eco_feel$upper.random),
  Lower90 = c(m.gen_eco_feel$TE - (1.645 * m.gen_eco_feel$seTE), m.gen_eco_feel$TE.random - (1.645 * m.gen_eco_feel$seTE.random)),
  Upper90 = c(m.gen_eco_feel$TE + (1.645 * m.gen_eco_feel$seTE), m.gen_eco_feel$TE.random + (1.645 * m.gen_eco_feel$seTE.random)),
  Narrative = "Economic Narrative",
  Outcome = "Feelings"
)

# Extract data for humanitarian treatment, deportation outcome
total_observations <- sum(m.gen_hum$n.e)
forest_data_hum <- data.frame( 
  Study = c(m.gen_hum$studlab, "Pooled Effect"), 
  Coefficient = c(m.gen_hum$TE, m.gen_hum$TE.random), 
  Observations = c(m.gen_hum$n.e, total_observations),
  SE = c(m.gen_hum$seTE, m.gen_hum$seTE.random),
  Lower95 = c(m.gen_hum$lower, m.gen_hum$lower.random),
  Upper95 = c(m.gen_hum$upper, m.gen_hum$upper.random),
  Lower90 = c(m.gen_hum$TE - (1.645 * m.gen_hum$seTE), m.gen_hum$TE.random - (1.645 * m.gen_hum$seTE.random)),
  Upper90 = c(m.gen_hum$TE + (1.645 * m.gen_hum$seTE), m.gen_hum$TE.random + (1.645 * m.gen_hum$seTE.random)),
  Narrative = "Humanitarian Narrative",
  Outcome = "Deportation"
)

# Extract data for humanitarian treatment, feelings outcome
total_observations <- sum(m.gen_hum_feel$n.e)
forest_data_hum_feel <- data.frame( 
  Study = c(m.gen_hum_feel$studlab, "Pooled Effect"),
  Coefficient = c(m.gen_hum_feel$TE, m.gen_hum_feel$TE.random),
  Observations = c(m.gen_hum_feel$n.e, total_observations),
  SE = c(m.gen_hum_feel$seTE, m.gen_hum_feel$seTE.random),
  Lower95 = c(m.gen_hum_feel$lower, m.gen_hum_feel$lower.random),
  Upper95 = c(m.gen_hum_feel$upper, m.gen_hum_feel$upper.random),
  Lower90 = c(m.gen_hum_feel$TE - (1.645 * m.gen_hum_feel$seTE), m.gen_hum_feel$TE.random - (1.645 * m.gen_hum_feel$seTE.random)),
  Upper90 = c(m.gen_hum_feel$TE + (1.645 * m.gen_hum_feel$seTE), m.gen_hum_feel$TE.random + (1.645 * m.gen_hum_feel$seTE.random)),
  Narrative = "Humanitarian Narrative",
  Outcome = "Feelings"
)


### Deportation plot ###

forest_data_combined <- rbind(forest_data_eco, forest_data_hum)

desired_order <- c("Pooled Effect", "Study 3", "Study 2", "Study 1")

forest_data_combined <- forest_data_combined %>%
  mutate(Study = factor(Study, levels = desired_order))

coefficient_plots <- ggplot(forest_data_combined, aes(x = Coefficient, y = Study)) +
  geom_point(aes(shape = "95% CI", color = Study), size = 3) +
  geom_errorbarh(aes(xmin = Lower95, xmax = Upper95, color = Study), height = 0, linewidth = 0.9) +
  geom_errorbarh(aes(xmin = Lower90, xmax = Upper90, color = Study), height = 0, linewidth = 1.9) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "", x = "", y = "", caption = "ATE Open Migration") +
  theme_minimal() +
  theme(
    axis.title.x = element_blank(), 
    plot.caption = element_text(size = 11, hjust = 0.5),
    strip.text = element_text(size = 12),
    legend.position = "none"
  ) +
  scale_x_continuous(limits = c(-0.1, 0.3)) +
  scale_color_manual(values = c("Study 1" = "grey50", "Study 2" = "grey50", "Study 3" = "grey50", "Pooled Effect" = "red")) +
  facet_grid(~Narrative)


### Feelings plot ###

forest_data_combined_feel <- rbind(forest_data_eco_feel, forest_data_hum_feel)

forest_data_combined_feel <- forest_data_combined_feel %>%
  mutate(Study = factor(Study, levels = desired_order))

coefficient_plots_feel <- ggplot(forest_data_combined_feel, aes(x = Coefficient, y = Study)) +
  geom_point(aes(shape = "95% CI", color = Study), size = 3) +
  geom_errorbarh(aes(xmin = Lower95, xmax = Upper95, color = Study), height = 0, linewidth = 0.9) +
  geom_errorbarh(aes(xmin = Lower90, xmax = Upper90, color = Study), height = 0, linewidth = 1.9) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(title = "", x = "", y = "", caption = "ATE Positive Affect") +
  theme_minimal() +
  theme(
    axis.title.x = element_blank(), 
    plot.caption = element_text(size = 11, hjust = 0.5),
    strip.text = element_text(size = 12),
    legend.position = "none"
  ) +
  scale_x_continuous(limits = c(-0.1, 0.3)) +
  scale_color_manual(values = c("Study 1" = "grey50", "Study 2" = "grey50", "Study 3" = "grey50", "Pooled Effect" = "red")) +
  facet_grid(~Narrative)


### Combined plot ###

combined_plot <- coefficient_plots / coefficient_plots_feel

combined_plot

ggsave(filename = paste0(plot_path, "figure_01.png"), plot = combined_plot, width = 7, height = 6, dpi = 600)